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Abstract. 

I overview the current understanding of the evolution of massive 
black hole (MBH) binaries in the center of the host stellar system. One 
jy^ . of the main questions is whether the stellar dynamical effect can make 

in 



the MBH binary hard enough that they can merge through gravitational 
wave radiation. So far, theories and numerical simulations suggested 
otherwize, since the hardening timescale becomes very long once "loss 
cone" is depleted. I'll present the result of recent simulations on this 
hardening timescale, and discuss its implication on the formation history 
of massive black holes. 
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There are now plenty of evidences that many, if not most, of giant ellipticals con- 
tain massive central black holes (Magorrian et al., 1998). Also, it has been sug- 
gested that the black hole mass Mbh shows tight correlation with the spheroidal 
mass and the central velocity dispersion (Gebhardt et al., 2000; Ferrarese & 
Merritt, 2000). The most straightforward explanation of such correlation is that 
massive galaxies are formed by merging of less massive galaxies, and that the 
central black hole also grows through merging (Kauffmann &: Haehnelt, 2000). 
This merging scenario naturally explains the observed correlations. 

This merging scenario has additional important characteristic that it nicely 
explains the observed structure of the central region of massive elliptical galax- 
ies. High-resolution observations by HST revealed that the "cores" of the giant 
elliptical galaxies are not really cores with flat volume density, but very shallow 
cusps expressed as p oc r~ a , where p is the volume luminosity density and r is 
the distance from the center, with power index a = 0.5 ~ 1 (Gebhardt et al., 
1996; Byun et al., 1996). Makino & Ebisuzaki (1996) performed the simulation 
of repeated mergers of galaxies with central black holes. They found that the 
structure of the merger converges to one profile. The merger product has a cen- 
tral cusp around the black hole with slope approximately —0.5, and the total 
mass of the stars in the cusp region is around the mass of the black hole binary. 
Nakano & Makino (1999a, b) showed, by combination of simple iV-body simula- 
tion and an analytic argument, that this shallow cusp is explained by the fact 
that distribution function of stars has a lower cutoff energy Eq. Their analytic 
argument naturally explains the observed correlation between the cusp radius 
and the black hole mass. 
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One remaining problem is what will happen to the binary black hole. Begel- 
man et al. (1980, hereafter BBR) predicted that the hardening of the BH binary- 
would be halted once it ejected all the stars in its "loss cone". Early numerical 
simulations (Makino et al., 1993; Mikkola & Valtonen, 1992) could not cover 
the range of the number of particles (relaxation time) wide enough to see the 
change in the timescale. Recent studies (Makino, 1997; Quinlan & Hernquist, 
1997; Chatterjee et al., 2003; Milosavljevic & Merritt, 2001) reported the re- 
sults, which are not only mutually inconsistent but also inconsistent with the 
theoretical prediction. BBR predicted that after the stars in the loss cone is de- 
pleted, the hardening timescale would be the relaxation timescale of the parent 
galaxy, since in the relaxation timescale stars would diffuse into the loss cone. 
The results of numerical simulation ranges from no dependence on the relaxation 
timescale (Milosavljevic & Merritt, 2001), to some dependence with upper limit 
in the timescale (Quinlan & Hernquist, 1997; Chatterjee et al., 2003), and finally 
to weaker-than-linear dependence (Makino, 1997). 

Thus, it has been an open question whether or not the loss-cone depletion 
occurs. Some even argued that numerical simulations would not help in deter- 
mining the fate of a BH binary. Here, we present the result of our resent iV-body 
simulations with up to 1M particles, where we saw a clear sign of the loss-cone 
depletion. 

2. Numerical method and Initial Models 

The details of the model is given elsewhere(Makino & Funato, 2003). We 
performed iV-body simulations using direct summation method on GRAPE- 
6(Makino et al., 2003). Gravitational interaction between field particles is soft- 
ened, while that involves BH particles is not. 

The initial galaxy model is a King model with nondimensional central po- 
tential of Wq = 7. We use the Heggie unit, where the mass M and the virial 
radius R v of the initial galaxy model and gravitational constant G are all unity. 
The mass of black hole particles is Mbh = 0.01. They are initially placed at 
(±0.5, 0, 0) with velocity (0, ±0.1, 0). 

The largest calculation (1 million particles and up to t = 400) took about 
one month, on a single-host, 4-processor-board GRAPE-6 system with the peak 
speed of 4 Tflops. 

3. Results 

Figure 1 shows the time evolution of the specific binding energy (per reduced 
mass) -E5 of the black hole binary. It is clear that the hardening timescale 
becomes longer as we increase N (and thereby the relaxation time). 

To quantitatively evaluate the dependence of the hardening rate on the 
number of particles, we calculated the hardening rate (3 (see figure 1). When 
Ef> = —1, the hardening rate is almost independent of N. However, as the binary 
becomes harder, (3 decreases, and the decrease is larger for larger N. Thus, the 
hardening rate (3 for large values of \Eb\ shows strong dependence on the number 
of particles N. For E^ = —7, the hardening rate (3 is almost proportional to 
1/N. In other words, we obtained the result which is consistent with the theory. 
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Figure 1. (Left) the evolution of the specific binding energy Ef, of the 
binary black hole. Curves are result with N = WK, 20K, 50K, 100K, 
200K, 500K and 1M (left to right). For all calculations, the softening 
length is e = 0.01. (Right) the hardening rate (3 = |A22;,/Ai| plotted as 
a function of the number of particles N. Dotted, dash-dotted, dashed 
and solid curves denote the values measured from the time at which 
| .Eft | reached 1, 3, 5 and 7, respectively. The end time is the time at 
which \Ef,\ reached initial value +0.5. 

The hardening timescale is proportional to the relaxation timescale of the parent 
galaxy. 

4. Discussions 

Our present result is consistent with the theory of loss-cone depletion and its 
refilling in thermal timescale, while previous numerical results are not. What 
caused this difference? 

Makino (1997) obtained the dependence much weaker than linear from sim- 
ulations with N = 2, 048 to 266,144. If we compare his figure 1 and our figure 1, 
the reason is obvious. The value of the binding energy at which the dependence 
is measured was rather small. In other words, Makino (1997) did not cover long 
enough time. 

Quinlan & Hernquist (1997) performed simulations very similar to the ones 
presented here. They found that the hardening rate was practically the same 
for N = 10 5 and N = 2 x 10 5 . We do not really understand why they obtained 
this result. Thy used the SCF method(Hernquist & Ostriker, 1992) and they 
varied the masses of particles depending on their initial angular momenta. We 
suspect this combination might have complicated the dependence of the hard- 
ening timescale on the number of particles. Chatterjee et al. (2003) might have 
similar problem. 

Milosavljevic & Merritt (2001) performed three runs with 8K, 16K and 
32K stars, and found that the hardening rate was independent of the number 
of particles. As suggested in their paper, this result is simply because the loss 
cone was not depleted in their simulation, primarily because N was too small. 
So their result is not inconsistent with ours. 
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As first suggested by BBR and confirmed by a number of foliowup works, 
if the hardening timescale of the black hole binary is proportional to the relax- 
ation timescale of the parent galaxy, the evolution timescale of a typical binary 
black hole in an elliptical galaxy exceeds the Hubble time by many orders of 
magnitude. In other words, the binary is unlikely to merge through encounters 
with field stars and gravitational wave radiation. 

Our result strongly suggests that the hardening timescale is indeed deter- 
mined by the relaxation timescale, for large enough N and after the binary 
becomes sufficiently hard. Thus, our result imply that gravitational interaction 
with field stars is insufficient to let the binary merge. 

If the binary black hole has long lifetime, it is quite natural to assume that 
some of the host galaxies which contain binary black holes would further merge 
with another galaxy with a central massive black hole or a binary. Simple esti- 
mate assuming the thermal distribution of the eccentricity (Makino & Ebisuzaki, 
1994) suggests that, during repeated three body interactions, the eccentricity of 
the binary can reach a very high value, resulting in quick merging through grav- 
itational wave radiation. Thus, strong triple interaction of three supermassive 
black holes may be common. 
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